#############################################################################
# R Script for Analysis of the AAA Motels Survey                            #
#############################################################################
#
setwd("d:/data") #Set the path
#
library(survey)
#
#
# Read in data and restrict to those with non-missing weights.
# Note that these are DRG replicate weights for G=10 dependent
# random groups, so weights are original weights times 10.
#
aaa<-read.table("AAA_motels_example.txt",header=TRUE)
aaa<-subset(aaa,WT!="NA")
#
# Construct intercept ("motel") and indicators for how frequently people ask 
# motel to make reservations for them.
# 
motel=rep(1,length(aaa$WT))
freqly=rep(0,length(aaa$WT))
never=freqly
rarely=freqly
freqly[aaa$RESERVAT==1]<-1
rarely[aaa$RESERVAT==2]<-1
never[aaa$RESERVAT==3]<-1
numer<-rarely+never
denom<-numer+freqly
rarely
#
# Add new variables to the data frame.
aaa<-data.frame(aaa,rarely,freqly,never,numer,denom,motel)
#
# Now compute standard errors using the DRG method, in two different ways: 
# directly, and using the R survey package.
# Directly:
G<-10
f.hat<-rep(0,G)
r.hat<-rep(0,G)
n.hat<-rep(0,G)
for(g in 1:G){
	f.hat[g]<-sum(freqly[aaa$RG==g]*aaa$WT[aaa$RG==g])
	r.hat[g]<-sum(rarely[aaa$RG==g]*aaa$WT[aaa$RG==g])
	n.hat[g]<-sum(never[aaa$RG==g]*aaa$WT[aaa$RG==g])
}
# Compare the following output to Table 2.3.3 in the Introduction to Variance
# Estimation.
cbind(1:G,f.hat,r.hat,n.hat)
#
#
numer.hat<-r.hat+n.hat
denom.hat<-numer.hat+f.hat
# Now compute variance and standard error for the 
# nonlinear statistic, (numerator)/(denominator):
mean(numer.hat/denom.hat)
var(numer.hat/denom.hat)/10
sqrt(var(numer.hat/denom.hat)/10)
#
# Second method, using the R survey package:
# Specify the survey design.
aaa.design<-svydesign(~ID,weights=~WT,data=aaa)
#
# Note that DRG weights are 10 times the survey weights, so divide by 10
# to get "Parent Sample" estimate. Ignore SE=standard error.
svytotal(~freqly+rarely+never,aaa.design)
#
# Table 2.3.3 (ignore standard errors here).
svyby(~freqly+rarely+never,by=~RG,aaa.design,svytotal)
#
# Ratio of �Rarely� plus �Never� to �Frequently� plus �Rarely� plus �Never�
#
RG.ratios<-svyby(~numer,denominator=~denom,by=~RG,aaa.design,svyratio)
# Compute the mean and DRG standard error estimate. (Note division by 10).
mean(unlist(RG.ratios[,2]))
var(unlist(RG.ratios[,2]))/10 # 
sqrt(var(unlist(RG.ratios[,2]))/10)
#
# Other exercises: go back and compute DRG estimates and standard errors 
# for the number of motels that frequently, rarely, or never 
# make reservations.
#
q() #Exit R